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ABSTRACT. We have used molecular dynamics to calculate the thermal conductivity of symmetric 
and asymmetric graphene nanoribbons (GNRs) of several nanometers in size (up to ~4-nm wide and 
~10-nm long). For symmetric nanoribbons, the calculated thermal conductivity (e.g. -2000 W/m-K 
@400K for a 1.5 nm x 5.7 nm zigzag GNR) is on the similar order of magnitude of the experimentally 
measured value for graphene. We have investigated the effects of edge chirality and found that 
nanoribbons with zigzag edges have appreciably larger thermal conductivity than nanoribbons with 
armchair edges. For asymmetric nanoribbons, we have found significant thermal rectification. Among 
various triangularly-shaped GNRs we investigated, the GNR with armchair bottom edge and a vertex 
angle of 30° gives the maximal thermal rectification. We also studied the effect of defects and found 
that vacancies and edge roughness in the nanoribbons can significantly decrease the thermal 
conductivity. However, substantial thermal rectification is observed even in the presence of edge 
roughness. 
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Graphene, a monolayer of graphite and two dimensional honeycomb lattice of sp bonded carbon, is 
the building block of most carbon based nanomaterial, such as carbon nanotubes (CNTs) and 
buckyballs. Much attention has been given to the exceptional and unique electronic properties of 
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graphene " uncovered in the past few years. Graphene nanoribbons (GNRs), narrow strips of graphene 
with few or few tens of nm in width, are particularly interesting and have been considered as important 
elements in future carbon-based nanoelectronics. For example, it has been shown that many of the 
electronic properties of GNRs may be tuned by its width or edge structures. 6 " 9 

In addition to its electronic properties, the thermal properties of graphene are also of fundamental and 
practical importance. Recent experiments 10 ' 11 have demonstrated that graphene has a superior thermal 
conductivity, likely underlying the high thermal conductivity known in CNTs 12 ' 13 and graphite 14 (ab- 
plane). This opens numerous possibilities of using graphene nanostructures for nanoscale thermal 
management. In this work, we have used molecular dynamics simulation to study thermal transport in 
few-nanometer-sized GNRs (typically 1-4 nm wide and 6-10 nm long). We have observed that the 
thermal conductivity of GNRs depends on the edge-chirality and can be affected by defects. We have 
also observed thermal rectification (TR) in asymmetric GNRs, where the thermal conductivity in one 
direction is significantly different from that in the opposite direction. Although TR has been 
experimentally observed in asymmetrically mass loaded nanotubes 15 and theoretically predicted in 
several other carbon nanostructures such as carbon nanohorns 16 and carbon nanotube intramolecular 
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junctions , it has not been studied in any graphene systems to our knowledge. TR have potential 
applications in nanoscale thermal management such as on-chip cooling and energy conversion by 
controlling the heat transport, and is also fundamental in several recently proposed novel schemes of 
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"thermal circuits" or information processing using phonons. " Our investigations of the thermal 
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conductivity and controlling heat flow in graphene nanostructures can be important for the development 
of energy-efficient nanoelectronics based on graphene. 
In our work, we have used classical molecular dynamics (MD) simulation based on the Brenner 

22 

potential of carbon-carbon interaction. We place atoms at the two ends of a GNR in the Nose-Hoover 
thermostats 23 ' 24 with temperatures T L (left end) and T R (right end) respectively (the temperature 
difference is denoted as A7"), and calculate the resulting heat current J . The thermal conductivity k is 
calculated from the well-known Fourier's law rc = Jd lijSTwh) where d , w and h (=0.335nm) are the 
length, width and thickness of the GNR respectively. The equations of motion for atoms in either the 
left or right Nose-Hoover thermostat are: 



d _ d _ 1 

dt ' dt t' 
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where the subscript i runs over all the atoms in the thermostat, p t is the momentum of the z'-th atom, F i 
is the total force acting on the z'-th atom, y and z are the dynamic parameter and relaxation time of the 
thermostat, T{t) is the instant temperature of the thermostat at time t , T (= T L or T R ) is the set 
temperature of the thermostat, N is the number of atoms in the thermostat, k B is the Boltzmann 
constant and m is the mass of the carbon atom. These equations of motion are integrated by a 3rd-order 
prediction-correction method. The time step is 0.5 fs, and the simulation runs for 10 7 time steps giving a 
total MD time of 5 ns. The relaxation time r is set to be 1 ps (see Additional Note 3). Typically, T{t) 
can stabilize around the set value T a after 0.5 ns. Time averaging of the temperature and heat current is 
performed from 2.5 ns to 5 ns. The heat current injected to the thermostat is given by 

J = Z,(- )pi 1 m = -wfc B T(t) ■ 
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Figure 1 . Structures of GNRs in this study: symmetric rectangle (a) and asymmetric triangle (b) and 

trapezoid (c). 

In Figure 1 , the atoms denoted with (black) squares at the ends are fixed to avoid the spurious global 
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rotation of the GNRs in the simulation. We have also performed the simulations with free and periodic 
boundary conditions, and found that the conclusions do not change qualitatively. The atoms denoted 
with triangles are placed in the Nose-Hoover thermostats (obeying equation (1)) at temperature T = T L 

(for the left-pointing triangles) and T = T R (for the right-pointing triangles) respectively. The atoms 

denoted with circles obey the Newton's law of motion — p. = F t , The average temperature of the two 

dt 

thermostats is T = (T L +T R )/2 . We define T L =(l±a)T and T R =(l + a)T ( a > ) so that the 
temperature difference is A7" = \T L — T R = 2aT . The parameter a - 10% is used in all of the following 

calculations (see Additional Note 3). At steady state, there is no energy accumulation in GNRs, so in 
principle J L +J R = . However, due to the nature of the discretized MD simulation, J L + J R slightly 
fluctuates around zero. Therefore, we define the heat current as J = (j L -J R )/2 and the error of the 
heat current as AJ — \J L + J R \ 1 2 (corresponding to the error in the thermal conductivity 
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Ak = AJd /(ATwh)). The temperature calculated from the classical MD (T MD ) has been corrected by 
taking into account the quantum effects of phonon occupation, using the scheme 

2T 3 ct d it x 2 

T MD = — y\ ~ dx , where T D (322 K) and T are the Debye temperature and corrected 
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temperature respectively (see Additional Note 3). 

First, we calculate the thermal conductivity of symmetric GNRs of rectangular shapes. For a 5.7-nm 
long and 1.5-nm wide rectangular GNR (with zigzag long edge), the thermal conductivity is found to be 
around 2100 W/m-K at 400 K (Figure 2a), on the same order of magnitude with the experimental 
measured value (-3000-5000 W/m-K) of graphene. 10 ' 11 On the other hand, the calculated thermal 
conductivity is nearly doubled after doubling the length of GNRs (see Additional Note 3). This suggests 
that our calculated thermal conductivity is limited by the finite length of GNRs and not corresponding to 
the value for graphene of macroscopic size. This is consistent with the phonon mean free path (MFP) in 
graphene extracted from the experiment (775 nm) 11 being much larger than the length (up to -10 nm) of 
the GNRs simulated in this study. We have also found that the calculated thermal conductivity remains 
nearly the same after doubling the width of GNRs (with the length unchanged). As we will see for all of 
the GNRs in this study, the thermal conductivity always monotonically increases with temperature (7), 
in the range we studied (100-400 K). Similar behavior is predicted in a recent theory on the thermal 
conductivity of small graphene flakes. 

The effect of edge chirality on the thermal conductivity in rectangular GNRs is also investigated. The 
chirality of GNRs (see the right inset of Figure 2b for the definition of chirality angle) is defined 
according to the edge parallel to the long direction of the GNR (different from the convention for 
CNTs). In Figure 2a, thermal conductivity is plotted as a function of temperature for both zigzag 
(dashed line) and armchair (solid line) edged GNRs. We show that the thermal conductivity of zigzag 
GNR (dashed line) is 20-50% larger than that of the armchair GNR (solid line). Recent work by Jiang et 
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al. obtained qualitatively similar results as ours using an approach based on ballistic phonon transport. 
Similar effects have also been studied for CNTs, but there have been no generally accepted agreement 
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on the preferred chirality of CNT for heat conduction. ' CNTs are periodic in the azimuthal direction 
and they can be considered as GNRs with infinite width. All phonon modes can propagate in CNTs. 
However, in GNRs with finite width the phonon modes propagating along the heat flow direction 
dominate. We speculate that the difference between the thermal conductivity of armchair and zigzag 
GNRs is mainly due to the different phonon scattering rates at the armchair and zigzag edges and the 
effect of the finite size of GNRs. Figure 2b shows the thermal conductivity as a function of the chirality 
angle at the temperature of 400 K. Two peaks of the thermal conductivity at chirality angles of 0° 
(corresponding to armchair edge) and 30° (zigzag edge) can be clearly identified. The peak at 30° is 
higher than that at 0° as seen in Figure 2a. We only need to study the chirality angle from 0° to 60° due 
to the 6-fold rotational symmetry of graphene. The GNRs with chirality angle not equal to integer 
multiples of 30° have irregular edges. Phonons can be strongly scattered by these irregular edges, likely 
resulting in the relatively low thermal conductivity at these angles (thus the two peaks at 0° and 30°) as 
seen in Figure 2b. 




Chirality Degree 

Figure 2. (a) Thermal conductivity of armchair and zigzag GNRs. (b) Thermal conductivity as a 
function of chirality angle (defined in right inset). The left two insets are two typical GNRs with edge 

chirality of zigzag (top) and armchair (bottom). 
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We have studied thermal conductivity of the asymmetric GNRs of triangular shapes (Figure lb) and 
found significant thermal rectification. It is shown in Figure 3 a that the thermal conductivity from the 
narrower (N) to the wider (W) end ( k n ^ w ) of such a triangular GNR is less than that from the wider to 

the narrower end (k w ^ n ). The thermal rectification factor, defined as rj = (rc w ^ N -k n ^ w )I k n ^ w is as 

large as 120% (at 7M80 K), as shown in the right inset of Figure 3a. We have also found similar 
thermal rectification in trapezoid-shaped GNRs (see Additional Note 4), but with smaller rectification 
factor. In contrary, in the symmetric rectangle GNRs there is no thermal rectification, i.e., the thermal 
conductivity form the left (L) to the right (R) is the same as that from the right to the left (dashed lines 
in Figure 3b) within the MD uncertainty. In calculating the thermal conductivity of asymmetric GNRs, 
the width (w) is taken as the width at the middle of GNRs. The temperature dependence of the thermal 
conductivity of the rectangular or trapezoidal GNRs is found to be similar to that of the rectangular 
GNRs. The overall thermal conductivity of asymmetric GNRs is lower than that of comparably sized 
symmetric GNRs. 
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Figure 3. Thermal conductivity of (a) triangular asymmetric and (b) rectangular symmetric GNRs. 
The right inset of (a) shows the thermal rectification factor rj as a function of temperature. 
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Previous studies of thermal rectifications (TR) " have suggested that TR originates from the 
interplay between structural gradient and lattice nonlinearity. A similar mechanism can underlie the TR 
in our asymmetric GNRs. We have investigated several geometric variations of the triangular GNRs to 
search for the structure with largest thermal rectification. To make better comparison, the solid lines in 
Figure 4a and 4b are the calculated thermal conductivity of asymmetric triangular GNR shown in Figure 
3a. This particular structure has an armchair bottom edge and a 30° vertex angle (defined in Figure lb). 
Figure 4a shows the thermal rectification of the GNRs with different vertex angles but the same bottom 
edge length and chirality (armchair). The dashed lines (dotted lines) are for the GNRs with vertex angle 
of 45° (60°). Compared with the GNR with vertex angle of 30° (solid line), the GNR with vertex angle 
of 45° (60°) has lower (higher) thermal conductivity in both directions, but both GNRs with vertex 
angles of 45° and 60° exhibit less thermal rectification (see the inset of Figure 4a). For the GNR with 
vertex angle of 45°, its hypotenuse has irregular edge and the phonon scattering at this edge likely 
decreases both the thermal conductivity and the thermal rectification. In triangular GNRs with armchair 
bottom edge, only the GNRs with vertex angle of 30° and 60° do not have irregular hypotenuse edges. 
At large vertex angles, the triangular GNRs approach the symmetric rectangular GNRs, which has zero 
thermal rectification. Since both irregular edge scattering and large vertex angles can decrease the 
thermal rectification, we suggest that a vertex angle of 30° may be optimal for thermal rectification. We 
have demonstrated in Figure 2 that the symmetric zigzag GNRs have larger thermal conductivity than 
that of the armchair GNRs. This is also true for asymmetric triangular GNRs. In Figure 4b (4c), the 
GNR with armchair bottom edge and vertex angle of 30° (60°) (solid lines) has smaller thermal 
conductivity than the GNR with zigzag bottom edge and vertex angle of 30° (60°) (dashed lines), but 
the former has larger thermal rectification. Overall, among various triangular GNRs we investigated, the 
one with armchair bottom edge and vertex angle of 30° has the largest thermal rectification. 
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Figure 4. Thermal conductivity of triangular GNRs: (a) dependence of thermal conductivity on vertex 
angles for triangular GNRs with armchair bottom edges and effect of edge chirality of the bottom edge 
on thermal conductivity at vertex angle of (b) 30° and (c) 60°. The thermal rectification factors versus 
temperature for the various structures simulated are shown in insets. 



In reality, GNRs inevitably have defects. We have studied the effect on the thermal conductivity of 
GNRs due to two types of defects: circular vacancies and edge roughness. In Figure 5a, we show the 
thermal conductivity of GNR with single (dotted line) and double (dashed line) circular vacancies and 
rough edges (dash-dot line). Here a single circular vacancy is created by removing all six carbon atoms 
of a hexagon. Compared to the perfect rectangular GNR (solid line), the thermal conductivity decreases 
after introducing circular vacancies. The edge roughness of the symmetric GNR can also decrease its 
thermal conductivity. This is in qualitative agreement with recent theoretical prediction by Nika et al. 
that the thermal conductivity of graphene flakes depends on the edge roughness and defect 
concentration, especially for small flakes. 30 The effect of edge roughness in the triangular GNR is also 
studied. We remove six atoms at the bottom edge and the hypotenuse of the triangular GNR of Figure 
lb. The solid lines in Figure 5b represent the thermal conductivity of the perfect triangular GNR with 
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armchair bottom edge and vertex angle of 30°. The edge roughness decreases both the thermal 
conductivity and the thermal rectification, which is nonetheless still as large as 80% at 7-180 K. 
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Figure 5. Effect of various defects on thermal conductivity of GNRs: (a) single and double vacancies 
and edge roughness for a symmetric rectangular GNR and (b) edge roughness for an asymmetric 
triangular GNR (the left inset shows the thermal rectification factor). 



In summary, we have used classical molecular dynamics to calculate the thermal conductivity of both 
symmetric and asymmetric GNRs of few nanometers in size. The calculated thermal conductivity of 
symmetric rectangular GNRs is on the same order of magnitude as the value expected for graphene, but 
with differences likely caused by the finite sizes of GNRs. The thermal conductivity is shown to depend 
on edge chirality and the zigzag edge GNRs is shown to have larger thermal conductivity than that of 
the armchair edge GNRs. We have demonstrated the thermal rectification effect in asymmetric 
triangular and trapezoidal GNRs. The triangular GNR with vertex angle of 30° and armchair bottom 
edge is found to have the largest thermal rectification among all GNRs studied. The defects can reduce 
the thermal conductivity of GNRs as well as the thermal rectification in asymmetric GNRs. 
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ADDITIONAL NOTES: 

32 

1) After this work was published, we became aware of similar results on thermal rectification in 
asymmetric GNRs obtained by Yang et al. 

2) The Fortran code used in this simulation is downloaded from 
http://www.eng.fsu.edu/~dommelen/research/nano/brenner/. The bugs of this code have been fixed by 
Dr. Leon van Dommelen. We have introduced Nose-Hoover thermostats into the codes and set the two 
thermostats at different temperatures to simulate the non-equilibrium thermal transport. 

3) More details about the MD calculations and the size dependence of thermal conductivity can be 
found in Ref. 34. 

4) More details about thermal rectification in trapezoid-shaped GNRs can be found in Ref. 35. 
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